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1. INTRODUCTION 



The most important sources for laser-interferometric gravitational-wave detec- 
tors like LIGO or VIRGO are catastrophic events such as coalescence of a 
neutron-star binary. A binary system evolves quasi-stationarily and the wave 
pattern is regular up to the onset of coalescence. In the last three minutes 
of coalescing binaries, comparison of observed wave patterns with theoretical 
templates will hence make it possible to determine mass and spin of stars. In 
this phase, the post-Newtonian expansion will give good templates since v ro %/c 
is as small as 0.1. In the final phase of the coalescence, however, there is no 
such small parameters. Numerical relativity is one of methods for studying 
such nonlinear phase of coalescing binary neutron stars; we call it the last three 
milliseconds. In order to know the characteristics of the waves in the last three 
milliseconds, general relativistic (GR) calculations are required. However, co- 
alescence of binary stars is a completely non-axisymmetric, 3 dimensional (in 
space) event and 3D GR calculation requires great powers of computers. We 
should hence make progress step by step. First of all, we used a Newtonian 
hydrodynamics code including radiation reaction of gravitational waves and 
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then we considered post-Newtonian correction, which is called (1+2.5) post- 
Newtonian calculation, since a radiation reaction potential is expressed as 2.5 
post-Newtonian terms. These calculations give a lot of perspectives on gravi- 
tational wave patterns and dynamics of coalescing events. Finally we began to 
attack a 3D GR code for the final phase of coalescing binary neutron stars. 

In this lecture, we describe results of numerical simulations of coalescing 
binary neutron stars using Newtonian and post-Newtonian hydrodynamics code 
and then discuss recent development of our 3D GR code for the last three 
milliseconds. 



2. POST-NEWTONIAN SIMULATIONS 

2.1. (1+2.5) Post-Newtonian Hydrodynamics Equations 

The equations of motion including general relativistic effects up to order (v/c) 2 
is called the first post-Newtonian (1PN) approximation. However, gravita- 
tional damping effects start at order (v/c) 5 , second- and-a-half post-Newtonian 
(2.5PN) terms, and the formulation for adding the effects to the equation of 
motion depends on the gauge condition. For example, the radiation reaction 
is represented by adding to the standard Newtonian potential a "radiation- 
reaction" potential proportional to the fifth time derivatives of the quadrupole 
moments jjj. 

In numerical calculations, however, it is hard to calculate the fifth time 
derivatives with satisfactory accuracy. Blanchet, Damour & Schafer || pre- 
sented (1+2.5)PN hydrodynamics equations including radiation damping ef- 
fects, where we need up to the third time derivatives of the quadrupole mo- 
ments. Following them, we use the evolution equations given by 

dtp + djW) = 0, (1) 

d t ( P W i ) + d j {pW i V^) = 2^ + j? PN +i?°«, (2) 

d t {pe) + d 3 (pev 3 ) = -pdjV , (3) 

where p, e and p are the coordinate rest-mass density, the internal energy 
density and the pressure, respectively. In order to express a hard equation of 
state for a neutron star, we use a polytropic equation of motion 

P = (7 - l)pe (4) 

with 7 = 2. The quantity v % denotes the 3-velocity and Wi is the "momentum 
per unit rest-mass." In terms of the 3-velocity v l is given by 

„<= (i- + + gt^gm (5) 



where 



■ij 



[3] - STF <{ 2 / [pdiWj - 2pw i d j ij) + x l d j i>d k {pw k ) - px l d^ t \ dV \ , (6) 
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which is the third time derivative of the reduced quadrupole moment of the 
system with the addition of corrections of order 0(c~ 2 ) ||. Here STF is the 
operator which takes the symmetric, trace-free part of any two-index object 
A ij : 



STF \A tj } = -A ij 
1 1 2 



1 



A 3% 



1 



The forces are given by 

F pre SS= _^ 

F}™ = -p 



■diii> 2 



^2 W A A 3 



(7) 



(8) 



Here tp is the Newtonian potential, a, /3, 6, ip2 and Ai are 1PN quantities and 
ip r is the radiation-reaction potential (a 2.5PN quantity): 



a = 


(37- 


- 2 )V>- ^T™ 2 , 


= 


1 2 


+ 7£ - 3^, 


5 = 


3 2 


+ (37 - 2)e + ^, 


At = 




1 i i 


Ipr = 


2 G 





(9) 



where iu 2 = S^WiWj and ^ is the time derivative of ^ with the addition of 
corrections of order (3(c~ 2 ). To compute tp, ipt, ip2, R and Ui, we must solve 
seven Poisson equations at each time step: 



Alp = 4nGp, 

Aip t = -inGdj (pwj) , 

Atp 2 = 47tG/9(5, 

Ai? =4TrGQ [ §x i d j p, 

AUi = -A-kG (Apwi + i x^djipwj) 



The energy flux of the gravitational waves is given by 



dE__G_ [3] ± 
dt ~ 5c 5 W « di y: 



(10) 

(11) 
(12) 

(13) 
(14) 



(15) 
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where 

Iij = STF 1 2 J p [wiWj - tfdjip] dv\ , (16) 

which is the second time derivative of the reduced quadrupole moment with 
the addition of corrections of order 0(c~ 2 ). The standard quadrupole formula 
gives the amplitude hj^ of the gravitational waves in the transverse-traceless 
gauge § ^ 

= ~ [PimPjn ~ ^Pij Pmn J jTq j (1') 

where Q mn is the reduced mass quadrupole moment 



Q mn = STF | J px m x n dVj (18) 

and Pij — Sij — riiUj is the projection operator onto the plane transverse to 
the outgoing wave direction, m — x % jr. For consistency with post-Newtonian 
accuracy, we must take into account the relativistic corrections up to the rela- 
tive order (y/c) 2 . This requires the time derivatives of the mass octupole, mass 
2 4 -pole, current quadrupole and current octupole Q . However we neglect these 
corrections since the gravitational wave amplitude is evaluated with quantities 
at each time step and small truncation errors above 'Newtonian' accuracy do 
not accumulate. With this accuracy, we can use ly defined by Eq.(|l6|) in place 
of Qij . From Eq. ( |l7| ) , two polarizations are given by 

hx = (19) 
where Iy is the quantity in the orthonormal basis 

Igg = {ixx COS 2 <j) + lyy Sm 2 <j) + 2ixy sin (j) COS (j)) COS 2 9 

+I ZZ sin 2 9 — 2 (I xz cos + I yz sin (f>) sin 9 cos 9, 
1^ = I xx sin 2 </> + lyy cos 2 <p - 2I xy sin <fi cos <fi, (20) 



Ig^ — (lyy — I xx ) cos 9 sin <p cos <p + Ixy cos #(cos 2 <p — sin 
+I XZ sin 9 sin </> — I yz sin 9 cos (f> 



2 



2.2. Numerical Methods 

The equations are discretized by the finite difference method (FDM) with a uni- 
form Cartesian grid. The evolution equations are integrated using van Leer's 
scheme j| with second-order accuracy in space. In order to make the scheme 
stable, a monotonicity condition, the so-called TVD (Total Variation Diminish- 
ing) limiter , is imposed on this scheme Q . In order to achieve second-order 
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accuracy in time we adopt a two-step procedure as follows. To describe our 
code, we write Eqs.(Q)-(^) as 

d t Q + d j (Qv^) = f(Q), (21) 

where Q is a 5-dimensional vector defined by 

Q=(p,p Wl ,pe). (22) 

The elements of 5-vector f(Q) are given by the right hand sides of Eqs.([l])- 
(H|). In order to obtain Q n+1 = Q(t n + Ai"), first we calculate the quantity 
Q n+ i = Q(t n + At n /2) from Q n as 

A in 

Q»+l = Q» + F(Q») _ , (23) 

and then calculate Q n+1 from Q n and Q n+ i as 

Q n+1 = Q n + F(Q n+ ^)M n , (24) 

where F(Q) is the flux given by j(Q). 

In order to treat shock waves, we use a tensor artificial viscosity given by 

( pe 2 (d k v k )STF{2d lV i} ifd k v k <0, 

Pii = \ . (25) 

I otherwise, 

where £ is an appropriate number with units of length. The gas pressure p 
is replaced by Py = pSij + pij . The Poisson equations are solved using the 
MICCG (Modified Incomplete Cholcsky-decomposition and Conjugate Gradi- 
ent) method 0, ||, |], 0, which is fully vectorized using the hyperplane 
method proposed by Ushiro [|ll| . 

We have performed various tests for this code JlJ, |l3), which include; (1) free 
transportation of a dust cube of a homogeneous or Gaussian density distribu- 
tion, (2) a ID Riemann shock tube, (3) a point explosion in the air, (4) local 
conservation of specific angular momentum for an axially symmetric collapse 
and (5) collapse of a homogeneous dust ellipsoid. Our code passed these tests 
with sufficiently good accuracy. 



2.3. Newtonian Calculation 

First we performed an extensive series of numerical calculations using Newto- 
nian hydrodynamics including the radiation damping effects, where terms of 
order 0(c~ 2 ) in Eqs.(|)-(|) are neglected while terms of order 0(c 5 ) are kept 
to include the radiation damping effects. 
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2.3.1. Coalescence of Neutron Stars without Spin 

We first consider initial data comprising of two neutron stars rotating rigidly 
around the center of mass. The binary system is assumed to be in rotational 
equilibrium at the initial time. 

In order to obtain hydrostatic equilibrium models of close neutron star bi- 
naries, we set dt = 0, _Fy cac = and v % — w\ in the Newtonian version of 
Eqs.([l])-(||). Here the gravitational radiation damping effects are neglected 
since a binary system evolves quasi-stationarily up to the onset of coalescence. 
Assuming the pressure p is given by p = Kp 1 and the two stars are in a synchro- 
nized circular orbit around the center of mass (xq, yo, 0), that is, the velocity is 
given by v % — (~(y~yo)Q, (x-xo)fl, 0), we have the equation that equilibrium 
models should satisfy: 



i, + H--{(x-x ) 2 + (y-yo) 2 }n 2 



0, 



(26) 



where H is the enthalpy, 



H = 



7- 



1 



Kp 



.7-1 



(27) 



A solution of Eq. (|2(|) for each star is given by 
1 



<p + H--{(x- x ) 2 + (y- y ) 2 } V 2 = C< (for i =1, 2) 



(28) 



where C\ and C 2 arc different constants in general. Since the position of the 
center of mass can be set freely, we set yo = 0. Instead of setting xo, however, 
we fix the positions x\ and x s 2 where the surface of each star intersects the x-axis 
between the stars; p(x\ , 0, 0) = p(x\, 0, 0) = 0. In addition, we fix the center of 
each star x\ and x\ and the density there; p(x^ 0, 0) = pi and p(x2, 0, 0) = p^. 
A self-consistent solution is determined by an iterative method. First, setting 
the constants 7 and K in the equation of state, we give an initial guess of the 
density distribution, usually as two spherical polytropcs. Then we repeat the 
following steps until convergence: 



(1) The potential ip is calculated as the solution of Eq.(10) 

(2) Ci, C2, £0 an d f2 are determined from 





-m- 


\& 


~xo) 2 n 2 


= Ci, 


Vx 




\{*\ 


~x Q ) 2 n 2 


= Ci, 




-m- 




~xo) 2 n 2 


= c 2 , 








~xo) 2 n 2 


= C 2 ; 



(29) 
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or explicitly 



A 2 {{x\) 2 -(x\) 2 }-A 1 {{x' i ) 2 -x{l) 2 } 
2{A 2 (x c 1 -xl)-A 1 (x c 2 -x s 2 )} 



^ = 7 \ , (30) 

{xl ~ xDixl + xl - 2x ) V ' 

where ^ a = i/>(x a , 0, 0),H a = H(x a , 0, 0) and A l = iff - + Hf. 
(3) Using these values, we determine a new density distribution from 

H = -^-Kp^ 1 =d-ip+ U(x- x ) 2 + y 2 }il 2 . (31) 
7 — 1 2 

Here we use the constant C\ for x > x$ and C 2 for x < xq. 

To be more realistic, we consider an infalling velocity due to the gravitational 
radiation. Assuming that the two stars are point masses with separation I in 
a circular orbit, I decreases at a rate 

^ = _ 64m 1 m 2 (m 1 +m 2 ) ^ 

We therefore add the infalling velocity given by this equation to the initial 
stars. 

We typically use a 141 x 141 x 131 grid. We assume reflection symmetry 
with respect to the z = plane and consider the region of z > only, since the 
system becomes rather flat and therefore needs a finer grid in the z-direction. 
Calculations were performed on a HITAC S820/80 supercomputer at the Na- 
tional Laboratory for High Energy Physics (KEK). Each calculation requires 
approximately 900Mbytes of memory. A typical CPU time required is 100-200 
hours for up to 50,000-90,000 time steps, which corresponds to an event lasting 
4-8 milliseconds. 

Here we show the results for neutron stars of 1.5M@ + 1.5M Q (EQSP2). 
Figure [l] shows the evolution of the density and the velocity on the x-y plane, 
where, at the initial time, the radius of each star is 8.8km, the separation 
between stars £q is 27km, the angular velocity f2 is 4.1 x 10 3 sec _1 . In each 
figure only the central part is represented, while the computational grid covers 
[—47km, 47km] in the x and y directions and [0km, 49km] in the z direction. 
Figure || shows the emitted energy L and the central density p c as functions 
of time. The gravitational radiation takes the angular momentum away from 
the system; in fact, a/m = J t /{GM 2 /c 3 ) = 0.64 at first and a/m = 0.38 
finally. Since the greatest part of the matter is within the Schwarzschild radius 
r = GM t /c 2 as shown in Fig.|l|i, the final destiny of the system must be a 
slowly rotating black hole. The energy emitted in gravitational radiation is 
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Fig. 1. — Density and velocity on the x-y plane for EQSP2. The time in units of 
milliseconds is shown. Solid lines are drawn in steps of a tenth of the maximum 
density and the inner and the outer dashed lines indicate 19/20 and 1/20 of maximum 
density. Arrows indicate the velocity vectors of the matter. A fat line shows a circle 
of radius 2GM t /c 2 to show the size of a spherical black hole for comparison. 
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Fig. 2. — Luminosity (in units of erg/sec) and central density (in units of g/cm ) 
functions of time (in units of millisecond) for EQSP2. 
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Fig. 3. - 



Wave forms of h + and h x observed on the z-axis at lOMpc for EQSP2. 
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Fig. 4. — Density and velocity on the x-y plane for TD1. 



1.6 x 10 53 erg, which is 3% of the total rest mass. This means that the efficiency 
of gravitational radiation is 30 times larger than the results for axisymmetric 
simulations by Smarr ju) or Stark and Piran |u| . 

Figure ^ shows the wave forms h+ and h x observed along the z-axis, which 
are given by 

^+ {I xx lyy) ■> 
2 

hx = ~Ixyi (33) 
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Fig. 5. — Luminosity (in units of erg/sec) and central density (in units of g/cm 3 ) as 
functions of time (in units of millisecond) for TD1. 



where r is the distance between the event and the observer. If the event occurs 
at lOMpc from the earth, the maximum amplitude of h is 4.0 x 10~ 21 . 

Figures M and | is the same as Figs.[I] and |], but for the model TD1, where 
the masses of stars are 1.70M© and 1.28M , the radii are 7km and 8km, 
SI = 7.9 x lf^sec- 1 and J t = 5.1GM|/c. 

2.3.2. Coalescence of Neutron Stars with Spin 

In the previous subsection, we place two neutron stars in a rotational equilib- 
rium state with rigid rotation around the center of mass as the initial condition. 
However, in the absence of the viscosity, the circulation T of the system should 
be conserved: 

r = / (V x v) ■ dS = f v-d£ = const. (34) 

JS JdS 

If the stars are assumed to rotate rigidly around the center of mass at first, the 
circulation in the z = constant plane is 2QqS, where f2o and S are the initial 
angular velocity and the area of the cross section of the stars, respectively. 
While two stars approach each other, the area of the cross section of the stars 
does not change so much, but the angular velocity J! becomes much larger than 
fio as SI oc (separation) 3 / 2 . Then, the stars must have spin with the angular 
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Fig. 6. — The same as Fig.|lp. The arrows show the velocity relative to the center of 
each star, which shows the spin retrograde to the orbital motion. 



velocity — fl to conserve the circulation. This is illustrated in Fig.|6|. It is the 
same figure as Fig.^je, but the arrows show the velocity relative to the center 
of each star. 

To take into account this effect, we consider two neutron stars rotating 
around the center of mass but also having spin angular velocity retrograde 
to the orbital motion fjqj . In this case, we don't know how to determine an 
equilibrium model because this is a similar problem to determine the Dedekind 
configuration p7| . As for an axisymmetric rotating star, however, we can ob- 
tain the equilibrium configurations. We here consider axisymmetric rotation 
stars as the initial condition for each spinning star. The equilibrium configu- 
ration for a spinning star is determined by Eq. (p9|) with ipi and f2 being the 
gravitational potential of each star alone and the spin angular momentum, re- 
spectively, as well as xo — x\. Then we assume that the orbital velocity of 
each star is the Keplarian one, VIk = y/G(mi + m^)/^, where mi, m2 and 
£0 are mass of each star and the separation of each star, and that ft = —fix. 
This initial condition is consistent if the tidal force is much smaller than the 
self-gravity. 

Initially the orbital velocity vk is given by 

v x K =-y£l K , v v K =xQ K (35) 
and the spin velocity v s is given by 



= y^K 
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Fig. 7. — Density and velocity on the x-y plane for SPIN1. 



x -^\n K (x>o) 

i (36) 
x+jjn K 0<o). 

Since v\ + = in the inertial frame, the velocity vector points to only 
y-direction. Figure |^ shows the evolution of the density and the velocity on 
the x-y plane, where, at the initial time, the mass and the radius of each star 
are 1.40M Q and 9km, the separation between stars £q is 27km. The evolution 
sequence is almost the same as EQSP2 up to the onset of coalescence. After 
coalescence starts, however, the evolution of the system is quite different; there 
is no spiral arm in the outer region and the configuration becomes a nearly 
axially symmetric disk soon. This is because the velocity in the outer region 
is smaller on account of the spin retrograde to the orbital motion and the 
centrifugal force is weaker than EQSP2. The stars are gradually coalescing 
because the centrifugal force is enhanced in the inner region, and the double 
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core structure is kept for a long time. The radiation reaction makes the double 
cores merge at last. 



2.4. Post-Newtonian Calculation 

Now we show results of a (1+2.5) post-Newtonian calculation we performed 
in order to study general relativistic effects. We performed also a Newtonian 
calculation with the same initial data for comparison. The initial data includes 
two spherical neutron stars, namely two relativistic polytropes with 7 = 2 
whose surfaces just touch each other. The mass and the radius of each star is 
0.62Mq and 15km, respectively. The two stars are rotating around the center of 
mass with angular velocity f2 = 2.0 x 10 3 sec _1 . The total angular momentum, 
J t , of the system is 1.6GMq/c. 

Figure | shows the evolution of the density and the velocity on the x-y plane 
and Figp shows the emitted energy L and the central density p c as functions of 
time. The result of the post-Newtonian(PN) calculation is shown on the right 
of Fig.|| and as solid lines in Fig.||, while the result of the Newtonian(N) calcula- 
tion is on the left and as dashed lines. In PN, coalescence begins more rapidly 
since general relativity effectively increases the gravitational force. Then a 
strong shock appears which makes the matter reexpand, and a radial oscil- 
lation with large amplitude is excited. The coalescence-expansion oscillation 
lasts longer in the PN case than in the N case. As a result the total energy 
emitted in the gravitational radiation increases in comparison with the N case; 
3.4 x 10 50 erg in PN and 2.2 x 10 50 erg in N up to 8.6msec. Figure [l0| shows the 
wave form observed along the z-axis, where solid lines are for PN and dashed 
lines are for N. Although h is damped faster in N, the maximum amplitudes 
are almost the same. The appearance of the strong shock in the PN calculation 
has little effect on the wave form. 

The difference between PN and N is caused principally by the strong shock 
which appears in the initial stage. Thus if the initial data consists of two stars 
in rotational equilibrium and coalescence proceeds more slowly, the difference 
between PN and N may be less pronounced. However we do not consider the 
red-shift effect and excitation of the quasi-normal modes of the resulting black 
hole or neutron star, which may be essential in the gravitational radiation 



1 15 1 Eq|. To complete the program of studying a final destiny of a coalescing 
binary system and computing emission of gravitational waves in such a event, 
we should solve a full set of the Einstein equation. 



3. GENERAL RELATIVISTIC SIMULATIONS 



In this section, we discuss recent developments of our fully general relativistic 
code for numerical calculation of the last three milliseconds of coalescing binary 
neutron stars. 
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Newtonian Post-Newtonian 




Fig. 8. — Density on the x-y plane. The left and right figures are for the Newto- 
nian(N) and post-Newtonian(PN) calculations, respectively. Notations are the same 
as for Fig.M. 
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Fig. 9. — Luminosity (in units of erg/sec) and central density (in units of g/cm 3 ) as 
functions of time (in units of millisecond) . The solid and dashed lines are for PN and 
N, respectively. 
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Fig. 10. — Wave forms of h+ and h x observed on the z-axis at lOMpc. The solid and 
dashed lines are for PN and N, respectively. 
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3.1. Basic Equations 

3.1.1. Initial Value Equations 

In the (3+l)-formalism of the Einstein equation, the line element takes the 
form 

ds 2 = -a 2 dt 2 + 7y (dx l + (3 i dt)(dx j + [3 J dt), (37) 

where a, (J 1 and jij are the lapse function, the shift vector and the intrinsic 
metric of 3-space, respectively. The initial data should satisfy the constraint 
equations 

R + K 2 - K i3 K lj = 16np H , (38) 
Dj(K* i -5 s i K)=&irJi, (39) 

where R is the 3-dimensional Ricci scalar curvature, Kij is the extrinsic cur- 
vature, K is its trace and the symbol Di denotes the covariant derivative with 
respect to 7y. The quantities p H and Ji are, respectively, the energy density 
and the momentum density of the matter measured by the normal line ob- 
server. The relation of these quantities to the stress-energy tensor is given by 
Eq.(H|) below. 

Now we assume that on the Cauchy slice the trace K = and 7y is confor- 
mally flat, that is, 

m = ^iih (40) 

where 7y is the flat metric. Defining the conformal transformation as 

= (j) 2 K^j , Kj = tfKj K^ = ^°K^ 
p B = (j) 6 p H and Ji = cjPJi, 

Eq. can be expressed as 

DjK't = torJi, (42) 
where Di is the covariant derivative with respect to 7^ p9). The traceless 



TT 

and the longitudinal traceless part (LW)ij [j20[; 



extrinsic curvature can be decomposed with the transverse traceless part K %] 



K l3 = + {LW) ih (43) 

where 

{LW)ij = DiWj + D 3 W t - I Jij&Wt. (44) 
Assuming Kj } T — 0, Eq.(ft2|) is reduced to 

AW, + ^DiD j W = 8^J l; (45) 
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where A = D l Di. Equation ( J45| ) is the coupled elliptic equation but, by 
defining 

X = DiW\ (46) 

it can be decoupled into 

A X = 6ttA J 1 (47) 

and 

AW t = 8nJ t - 1 D iX . (48) 
The boundary conditions for \ and Wi are 

X = O (i) (49) 

and 

e ijk x j M k „ / 1 



where M k is a constant related to the angular momentum of the system. 
The conformal transformation of the scalar curvature is 

R = <p~ 4 R- 8<?T 5 A0, (51) 

where R is the scalar curvature with respect to 7^. Since R = for a confor- 
mally flat metric, Eq.pSh is reduced to 



A(j> = -2n<f>- 1 p B - \<f> 7 A '*'. (52) 



The boundary condition for (f> is 

^ =1 + ^ + °(^)' (53) 
where Mq is the gravitational mass of the system. 
3.1.2. Relativistic Hydrodynamics 

We assume the perfect fluid stress-energy tensor, which is given by 

= (p + pe+ p)u^u v + pg^, (54) 

where p, e and p are the proper mass density, the specific internal energy and 
the pressure, respectively, and is the four-velocity of the fluid. The energy 
density p H , the momentum density Ji and the stress tensor Sij of the matter 
measured by the normal line observer are, respectively, given by 



n u 



T\ivi Ji — hi^n and Sij = hi^hj T^ v , (55) 
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where n M is the unit timelike four-vector normal to the spacelike hypersurface 
and /i M „ is the projection tensor into the hypersurface defined by 

= g^ v +Tfyn„. (56) 

The relativistic hydrodynamics equations are obtained from the conservation of 
baryon number, V M (pit A1 ), and the energy-momentum conservation law, V„T M y . 
In order to obtain equations similar to the Newtonian hydrodynamics equa- 
tions, we define p N and uf as 

Pn = Viau°p and uf = (57) 

respectively, where 7 = det("fij). Then the equation for the conservation of 
baryon number takes the form 

d t p N +d e (p N V e )=0, (58) 

where 

V = ^ = (59) 
u P + Ph 

The equation for momentum conservation is 

dt(p N u^) + d e (p N u?V e ) = -^/adtp - ^7(p + P H )d<a 

\ dau + ViW e - (60) 

2 {P + Ph) 

The equation for internal energy conservation becomes 

d t (p N e) + d e (p N eV e ) = -V^du (VW) . (61) 

To complete hydrodynamics equations, we need an equation of state, 

P = P(£,P)- (62) 

Equations (|58|)-(|6l"|), whose structure is the same as Eqs.(|l])-(|3]), can be solved 
using the post-Newtonian hydrodynamics code described the previous section. 

3.1.3. Time Evolution of the Extrinsic Curvature 

The evolution equation for the extrinsic curvature takes the form 

dtKij - d e {K l3 (i l ) = (SK)™ in + (SK)ff + (SK% , (63) 

where 

(SK)^ = a {R l3 - 8tt [Sij + | 7 y (p H - S e e )] } - D^a (64) 
{SK)ff = a (KKij - 2K u K e 3 ) , (65) 

{SK% = K mi d J p m + K mj d^ m - K l3 d m p m , (66) 

and Rij is the 3-dimensional Ricci tensor. Equation ([33]) can be solved in the 
same way as the hydrodynamics equations with the velocity — f3 e instead of V e . 
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3.1.4. Time Slicing 

Defining the conformal factor (j> as 

0=(det( 74J )p, (67) 
we use the lapse function a given by 

4>-if , 



q = exp 



-2 </>-! 



(68) 



Shibata and Nakamura pjj called this time slicing as the conformal time slicing 
because a is determined by the conformal factor cf>. In the conformal time slic- 
ing, the space outside the central matter quickly approaches the Schwarzschild 
metric. Detailed discussion on the nature of the conformal time slicing is given 
in the reference [^l] . The equation for the conformal factor <f> is obtained from 
the Hamiltonian constraint equations, 

A0 = -4- (!6np H + KijK ij -K 2 - ^R) , (69) 
o 

where A is the Laplacian with respect to 

lij = <l>~\j- (70) 
Note that 7 y is not the flat-space metric for t ^ 0. 
3.1.5. Spatial Coordinates 

Here we define conformal transformation, different from Eq. , as 



ff w = ^-*tf w , A',' A','. K = K, 



-4/ 



ft ee 0- 4 ft and ee /3* 
and the equation for metric thereby becomes 



(71) 



dtlij = f + Djfa - pijDrfA - 2a (k^ - hyyK j 1 72 ) 

where Dj denotes the covariant derivative with respect to 7^ defined by Eq. (JTC 
Now we demand the divergence with respect to the flat-space metric of the right 
hand side of Eq. ((72]) to vanish, 

diA% = 0, (73) 
which is reduced to the equation for the shift vector f3 l , 

V 2 /3' + \di (dtft) = 
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2a K., 



K 



(74) 



0, 



where V 2 is the simple flat-space Laplacian and jij = jij — Sij. Equation ( [74]) 
is the similar to Eq.(45) while (3 l appears in the right hand side, too. 
If we demand 

D j A% = 0, (75) 

instead of Eq.(7^), we have the minimal distortion condition proposed by York 
and Smarr |22|. In this sense, we call Eq.(f73|) the pseudo- minimal distortion 
condition. 



3.2. Numerical Methods 

We use Cartesian (x,y,z) isotropic coordinates for relativistic simulation, too. 
The hydrodynamics equations, Eqs.(|58|)-(|6l|), and the evolution equations for 
the extrinsic curvature, Eq. (|63|) , and the 3- metric, Eq.(|7^), are solved through 
the same scheme as in the Newtonian simulation. The elliptic equations are 
solved through the MICCG method. The equation for the conformal factor in 
the initial hypersurface, Eq.(|5|), can be reduced to 

V 2 = 4tt Pi (0), (76) 

since the initial hypersurface is conformally flat (A = V 2 ). The source term 
pi{4>) is a non-linear function of </>. This equation is solved by an iterative 
method: The following iteration is repeated until convergence, 



(i+i) - 



= (V 2 )- 1 [47T Pl (>>)] (far J =1,2,.-.). (77) 



Equations (|69|) and ( |74| ) for t > are also non-linear or coupled elliptic equa- 
tions, which can be written as 

V 2 = 47rp 2 (</>), (78) 

V 2 X = 3^(^(/3)) (79) 

and 

V 2 /? =~d iX + ^pi(P), (80) 

where x = de/3 e . The cost of finding self-consistent solutions for Eqs.(J7q)- 



(80) at each time step is very expensive and hence we use the values at the 
previous time step for <j> and fS 1 in the source terms p 2 and p 3 . The accuracy of 
solutions by this method seems to be sufficient to examine the global behavior 
of the numerical schemes and the coordinate conditions. The method should 
be improved in a future version. 

Although any equation of state can be used, for the present test calculation, 
we use a simple equations state Eq.(|j) with 7 = 2 to express a hard equation 
of state for a neutron star. 
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3.3. Numerical Results 

Now we show results of test calculations. Here we use units of 

M = M e , L = — ^ and T= — ^. 81 

c c J 

In the following test calculations, we use a 80 x 80 x 80 grid covering [—40, 40] 
in each direction. 

3.3.1. Spherically Symmetric Dust Collapse 

In the conformal time slicing given by Eq. ( |68| ) , the time- variable parts of 7^ 
can be considered as gravitational radiation parts, since the space outside the 
central matter approaches the Schwarzschild metric quickly [^IJ . Thus the total 
energy of the gravitational waves is given by 

2 



E GW = I dt I dO^-(4 T ) 2 , (82) 



where Aj^ is the transverse-traceless part of the time derivative of 



7ij, 



^ = (^7.) TT (83) 



and the time derivative of Ty is given by Eq.(|72j). Now we define the "energy 
density of the gravitational waves" , although it cannot be defined locally, as 

■w-sj^'-sf^"^™* (84) 

In order to examine whether Aj^ shows the gravitational waves' degree of 
freedom, we performed calculation of collapse of a spherically symmetric dust 
star of mass M. Figure |ll| shows the mass density p N and the "energy density 
of the gravitational waves" r 2 p Gw . Although r 2 p QW has a non-zero value near 
the center, it does not propagate outward. Since the meaning of "transverse- 
traceless" is clear only in the asymptotically flat region, the results is consistent 
with the idea that Ajj represents the gravitational waves. The total energy of 
the gravitational waves E GW is ~ 10~ 7 Af at t = 30M observed at r = 20M. 
This value is essentially zero within numerical errors. 

3.3.2. Formation of a Rotating Black Hole 

Stark and Piran ]lq | performed axially symmetric simulations for formation 
of a black hole and directly computed their gravitational radiation emission. 
For comparison with their results, we set similar initial conditions, although 
our coordinate condition is different from theirs. We use a simple equation of 
state, Eq.(||), which is the same as Stark and Piran. We place a spherically 
symmetric star of 7 = 2 polytrope of mass IMq with internal energy necessary 
for Newtonian equilibrium. The rotation velocity with constant f2 is added to 
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Density 
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MAX 

rp G ,.,= 5.0 x id " 



Fig. 11. — Mass density p N (left) and "energy density of the gravitational waves" 
r 2 p GW (right) on the x-y, y-z and x-z planes for a spherically symmetric dust collapse. 
The time in units of Mq, and the maximum value of density in units of Mq 2 are 
shown. 



this system. We performed four simulations, ROT1, ROT2, ROT3 and ROT4, 
each of which is characterize by the initial value of a/M g , where a is the angular 
momentum per unit gravitational mass M g . The value of a/M g is 0.3, 0.5, 0.8 
and 1.0, respectively. 

Figure [j^ show the evolution of density distribution of ROT4. It represents 
that the centrifugal force prevents the star from collapse in the direction per- 
pendicular to the rotation axis and the star becomes an oblate spheroid. The 
"energy of gravitational waves" is shown in Figjl^. A pattern almost axially 
symmetric with a few peaks appears on the x-y plane. Deviation from axial 
symmetry is caused by coarseness of our Cartesian grid. The wave pattern on x- 
z and y-z planes looks like the quadrupole emission proportional to sin 9, which 
is also the case in Stark and Piran. The total energy radiated at t = 30Af Q 
observed at r = 20M Q are 5.7 x 1O~ 7 M , 3.8 x 1O _6 M 0) 1.5 x 10" 5 M Q and 
3.6 x 10~ 5 Af o for ROT1, ROT2, ROT3 and ROT4, respectively. Stark and 
Piran found that the total energy of gravitational waves is proportional to a 
for small a and it levels off at a ~ M g . The same dependency reappears in our 
calculation; the total energy is expressed approximately 

£cw-5xl0- 5 ^j . (85) 

We have to remark that the value of the coefficient in Stark and Piran is about 
20 times larger than ours. However, this is not disagreement, because Stark 
and Piran pursued their calculation up to t = 100M while we did only up 
to t = 30M. The energy we observed is thus just expected from the time 
dependence of the gravitational wave energy given by Stark and Piran [[15| . 
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ROT4 
Density 




Fig. 12. — Mass density p N on the x-y, y-z and x-z planes for a spherically symmetric 
collapse of a rotating star R0T4. Arrows indicate the velocity vector V . 



3.3.3. Coalescing Binary Neutron Stars 

Now we begin to attack the problem of coalescing binary neutron stars. We 
place as initial data two spherical neutron stars of mass M — 1.0M© and 
radius r = 6M with density distribution (p N ) of 7 = 2 polytrope at x = r 
and x = — r ; two stars just touch each other. We add a rigid rotation with 
angular velocity ft as well as an approaching velocity v a to this system such 
that 

[ -ily + v a if x < 0, 
= n Xl (87) 

setting v a = ^r . 

We performed three simulations, BI1, BI2 and BI3, with different values of 
il. The total angular momentum divided by the square of the gravitational 
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Fig. 13. — "Energy density of the gravitational waves" r 2 ' p GW for ROT4. 



mass are 0.71, 1.43 and 1.01 for BI1, BI2 and BI3, respectively. 

Figure [l4| shows the evolution of the density on the x-y, y-z and x-z planes 
of BI3. The final object will be a rotating black hole although we have not 
yet tried to determine the apparent horizon. The "gravitational wave energy" 
p GW shows an interesting feature of generation and propagation of the waves. 
A spiral pattern appears on the x-y plane while different patterns with peaks 
around z-axis appear on the x-z and y-z planes. This can be explained naively 
by the quadrupole wave pattern given by 

2 

t 2 Pgw = Zr- (47) 2 « cos2 + sm 2 9sm 2 (2fl(t - r) - 2ip)/4. (88) 

327T 

On the x-y plane, where 9 = 7r/2, p GW is constant along the spiral of r + tp/Q = 
constant, while near z-axis, where 8 w 0, p GW cx cos 2 9. 

At the numerical boundary, we have to impose the outgoing wave condition 
to 7ij and Kij. In the present calculation, however, x max , 2/max and z m ax are 
comparable to the wave length of the expected gravitational waves, ~ 40Mq. 
We tried various boundary conditions, but none works well. We therefore adopt 
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a simple extrapolation at the numerical boundary and suspend calculations 
when the gravitational waves reach the numerical boundary. We hope that the 
problem will be overcome if the numerical boundary is located farther, where 
the outgoing wave condition becomes appropriate. 

The total energy of the gravitational radiation amounts to 3 x 1O~ 3 M at 
t = 35M„ for BI3. 
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Fig. 14. — Density on the x-y, y-z and x-z planes for BI3. Arrows indicate the 
velocity vector V . 
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Fig. 15. — Energy density of the gravitational waves r 2 p GW on the x-y, y-z and x-z 
planes for BI3. 
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